!
! Copyright (C) 2001-2009 Quantum ESPRESSO group
! Copyright (C) 2009 Brian Kolb, Timo Thonhauser - Wake Forest University
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
MODULE kernel_table
  
  !!  This module is used to read in the kernel table file
  !!  "vdW_kernel_table" and store some of the important parameters.  The
  !!  top of the vdW_kernel_table file holds the number of q points, the
  !!  number of radial points (r) used in the kernel generation, the maximum
  !!  value of r used (r is the parameter in the kernel function d=q*r where
  !!  q is defined in DION equation 11), and the values of the q points
  !!  used.  These parameters are stored as public parameters for use in
  !!  various routines.  This routine also reads the tabulated values of the
  !!  Fourier transformed kernel function for each pair of q values (see
  !!  SOLER equations 3 and 8).  Since these kernel functions need to be
  !!  interpolated using splines, the second derivatives of the Fourier
  !!  transformed kernel functions (phi_alpha_beta) are also tabulated in
  !!  the vdW_kernel_table and are read in here.

  !!  This is done in a module because there are quite a few subroutines in
  !!  xc_vdW_DF.f90 that require knowledge of the number (identity) of q
  !!  points, the maximum value of the radius, and, of course, the tabulated
  !!  kernel function and its second derivatives (for spline interpolation).
  !!  Putting this routine in a module meas that those routines can just use
  !!  kernel_table rather than passing variables around all over the place.
  
  USE kinds,                  ONLY : dp
  USE io_files,               ONLY : pseudo_dir
  USE constants,              ONLY : pi
  use wrappers,               ONLY : md5_from_file
  implicit none

  private
  save
  
  !! Variables to be used by various routines in xc_vdW_DF.f90, declared
  !! public so they can be seen from outside
  ! --------------------------------------------------------------------------
  
  public :: Nqs, Nr_points, r_max, q_mesh, q_cut, q_min, dk
  public :: kernel, d2phi_dk2
  public :: initialize_kernel_table
  public :: vdw_table_name
  public :: vdw_kernel_md5_cksum
  integer :: Nqs, Nr_points                      !! The number of q points and radial points
  !                                                    !! used in generating the kernel phi(q1*r, q2*r)
  !                                                    !! (see DION 14-16 and SOLER 3)

  real(dp) :: r_max, q_cut, q_min, dk            !! The maximum value of r, the maximum and minimum
  !                                                    !! values of q and the k-space spacing of grid points.
  !                                                    !! Note that, during a vdW run, values of q0 found
  !                                                    !! larger than q_cut will be saturated (SOLER 5) to
  !                                                    !! q_cut

  real(dp), allocatable :: q_mesh(:)             !! The values of all the q points used

  real(dp), allocatable :: kernel(:,:,:)         !! A matrix holding the Fourier transformed kernel function
  !                                                    !! for each pair of q values.  The ordering is
  !                                                    !! kernel(k_point, q1_value, q2_value)

  real(dp), allocatable ::  d2phi_dk2(:,:,:)     !! A matrix holding the second derivatives of the above
  !                                                    !! kernel matrix at each of the q points.  Stored as  
  !                                                    !! d2phi_dk2(k_point, q1_value, q2_value)
  !
  character(len=256) :: vdw_table_name = ' '           !! If present from input use this name
  CHARACTER(LEN=30)  :: double_format = "(1p4e23.14)"
  CHARACTER(len=32)  :: vdw_kernel_md5_cksum = 'NOT SET'
  !
  INTEGER, EXTERNAL :: find_free_unit
  ! --------------------------------------------------------------------------

CONTAINS


  !! ###################################################################################
  !!                          |                           |
  !!                          |  INITIALIZE_KERNEL_TABLE  |
  !!                          |___________________________|

  !!  Subroutine that actually reads the kernel file and stores the parameters.  This routine
  !!  is called only once, at the start of a vdW run.  

  subroutine initialize_kernel_table(inlc)
    
    integer, INTENT(IN) :: inlc

    integer :: q1_i, q2_i                           !! Indexing variables

    integer :: kernel_file                          !! The unit number for the kernel file

    logical :: file_exists                          !! A variable to say whether
    !                                               !! needed file exists

    character(len=1000) :: kernel_file_name         !! The path to the kernel file.
    !                                               !! Although this name must be
    !                                               !! "vdW_kernel_table", this variable
    !                                               !! is used to hold the entire path
    !                                               !! since we check 3 places for it.
    
    
    !!write(*,*) "Reading kernel table ... "

     !! Get the unit number for the kernel file
    kernel_file = find_free_unit()
   

    !!
    !!

    if (TRIM(vdw_table_name)==' ') then
      
      if (inlc==3) then
        vdw_table_name='rVV10_kernel_table'
      else
        vdw_table_name='vdW_kernel_table'
      endif   
    
    endif
 
   if (allocated(kernel)) then

       return

    end if

    
    !! First we check the current directory for the vdW_kernel_table file
    !! If it is not found there it is looked for in the pseudopotential 
    !! directory.  If it's not there the default kernel file installed
    !! in the PW directory of the PWSCF source is tried.  If none of those
    !! exist the code crashes.

    kernel_file_name=vdw_table_name
    inquire(file=kernel_file_name, exist=file_exists)

    !! If the file is found in the current directory we use that one
    !! ------------------------------------------------------------------------------------------

    if (.not. file_exists) then

       !! No "vdW_kernel_table" file in the current directory.  Try the pseudopotential directory
       !! -----------------------------------------------------------------------------------------

       kernel_file_name = trim(pseudo_dir)//'/'//vdw_table_name
       inquire(file=kernel_file_name, exist=file_exists)

       if (.not. file_exists) call errore('read_kernel_table', &
                         TRIM(vdw_table_name)//' file could be found',1)

    end if

    !! Generates the md5 file
    CALL md5_from_file(kernel_file_name, vdw_kernel_md5_cksum)

    !! Open the file to read
    open(unit=kernel_file, file=kernel_file_name, status='old', form='formatted', action='read') 

    !! Read in the number of q points used for this kernel file, the
    !! number of r points, and the maximum value of the r point
    read(kernel_file, '(2i5)') Nqs, Nr_points
    read(kernel_file, double_format) r_max

    allocate( q_mesh(Nqs) )
    allocate( kernel(0:Nr_points,Nqs,Nqs), d2phi_dk2(0:Nr_points,Nqs,Nqs) )
    
    !! Read in the values of the q points used to generate this kernel
    read(kernel_file, double_format) q_mesh


    !! For each pair of q values, read in the function phi_q1_q2(k).
    !! That is, the fourier transformed kernel function assuming q1 and q2
    !! for all the values of r used.
    !! ----------------------------------------------------------------------------------------------
    
    do q1_i = 1, Nqs
       do q2_i = 1, q1_i

          read(kernel_file, double_format) kernel(0:Nr_points, q1_i, q2_i)
          kernel(0:Nr_points, q2_i, q1_i) = kernel(0:Nr_points, q1_i, q2_i)

       end do
    end do

    !! ----------------------------------------------------------------------------------------------



    !! Again, for each pair of q values (q1 and q2), read in the value
    !! of the second derivative of the above mentiond Fourier transformed
    !! kernel function phi_alpha_beta(k).  These are used for spline
    !! interpolation of the Fourier transformed kernel.
    !! -----------------------------------------------------------------------------------------------

    do q1_i = 1, Nqs
       do q2_i = 1, q1_i

          read(kernel_file, double_format) d2phi_dk2(0:Nr_points, q1_i, q2_i)
          d2phi_dk2(0:Nr_points, q2_i, q1_i) = d2phi_dk2(0:Nr_points, q1_i, q2_i)

       end do
    end do

    !! -----------------------------------------------------------------------------------------------

    close(kernel_file)

    !! Define a few more vaiables useful to some of the subroutines in xc_vdW_DF.f90
    !! ------------------------------------------------------------------------------------------------

    q_cut = q_mesh(Nqs)
    q_min = q_mesh(1)
    dk = 2.0D0*pi/r_max

    !! ------------------------------------------------------------------------------------------------

  end subroutine initialize_kernel_table

  !! #################################################################################################



end MODULE kernel_table
